TP2 geneticData Mining - PCA-regression in genetics

Lamyaa BOUZBIBA Alexandre POUPEAU Eloise JULIEN

1. geneticData

NAm2 = read.table("NAm2.txt", header = TRUE)
# unique is used to get the name of all population / ethnics without duplicates
names=unique(NAm2$Pop)
npop=length(names)
# coordinates for each pop
coord=unique(NAm2[, c("Pop","long","lat")])
colPalette=rep(c("black","red","cyan","orange","brown","blue","pink","purple","darkgreen"),3)
# type 16 = circle / type 15 = square / type 25 = triangle
pch=rep(c(16,15,25),each=9)
# the display works because for each population we have the latitude and the longitude
plot(coord[, c("long","lat")],pch=pch,col=colPalette,asp=1)
# asp allows to have the correct ratio between  axis  longitude and latitude
# Then the map is not deformed  
legend("bottomleft",legend=names,col=colPalette,lty=-1,pch=pch,cex=.75,ncol=2,lwd=2)

library(maps);map("world",add=T)

The description of the script is in the commentaries above.

2. Regression

# turn the matrix into a geneticData frame
geneticData = NAm2[,-c(1:7)]
geneticDataFrame = data.frame(geneticData)
# creation of the linear model
reg = lm(formula = long ~ . , data = geneticDataFrame)
# uncomment that following part to see the summary of the regression
# summary(reg)

All values are NA (Not Available) because we have too much predictors. The \(X\) matrix is cannot be reversed so we have to use the PCA method to counter that problem.

3.PCA Principal Component Analysis

a)

The main goal of the PCA is to reduce the number of predictors. In order to achieve that, we have to extract the most important components of the geneticData set. These new variables are a linear combinaison of the original ones. The number of principal components is smaller than the number of original variables. The aim is to maximize the variation of the geneticData set.

b)

# install.packages("factoextra")
library(factoextra)

Loading required package: ggplot2
Welcome! Related Books: `Practical Guide To Cluster Analysis in R` at https://goo.gl/13EFCZ

allGenes <- geneticDataFrame[, -1]
pcaNAm2 <- prcomp(allGenes, scale = FALSE)
fviz_eig(pcaNAm2)

vectProportionOfVariance <- pcaNAm2$sdev^2/sum(pcaNAm2$sdev^2)
# summary(pcaNAm2.pca)

It is better to use scale = FALSE by convention because we get a better representation of the percentage of explained values.

FALSE : 494 = linear combianison of the 5000 genes We notice that we get only 494 principal components instead of approximately 5000 which correspond to the number of genes. This can be explained because “” in the covariance matrix the eigenvalues from 495 to ~5000 are equals to zero. In reality, we only select the genes that best explain the model up to 500. If we had to really find the “best-explaining”, which means finding the most explicative genes from 1 to 5000 we would have to test all the combinaisons of 494 genes.

c)

caxes=c(1,2)
plot(pcaNAm2$x[,caxes],col="white")
for (i in 1:npop) {
  # print(names[i])
  lines(pcaNAm2$x[which(NAm2[,3]==names[i]),caxes],
  type="p",col=colPalette[i],pch=pch[i])
}

legend("topleft",legend=names,col=colPalette,lty=-
1,pch=pch,cex=.65,ncol=3,lwd=2)

caxes=c(5,6)
plot(pcaNAm2$x[,caxes],col="white")
for (i in 1:npop) {
  # print(names[i])
  lines(pcaNAm2$x[which(NAm2[,3]==names[i]),caxes],
  type="p",col=colPalette[i],pch=pch[i])
}

legend("topright",legend=names,col=colPalette,lty=-
1,pch=pch,cex=.65,ncol=3,lwd=2)

The Ache population is well identified by both components. The Surui can also be identified by both components but more especially by the PC1.

The Karitiana population is well identified by the PC5. The Pima population is well identified by the PC6.

d)

inertia = vectProportionOfVariance[1] + vectProportionOfVariance[2]
sprintf("The inertia is equal to %f.", inertia)

[1] "The inertia is equal to 0.035684."

According to what we have seen until now, we have two ideas.

The first idea is that we have seen that we can distinguish two population using two principal components. Therefore, as there are 27 different populations, we could deduce that we would need 14 pairs of principal components.

However, we can ask ourselves if we can really each time distinguish two different populations using two principal components. It might not be the case because we have just tested it with two pairs of principal components. Thus, we could think about an alternative way, which is adding the variance of the first principal components until the sum is equal to a certain level. This level could be 25% or 50%. This last one means using a lot of principal components.

numberOfPCs <- function(vectorProportionOfVariance, percentage) {
  tmp = 0.0
  i = 1
  while (tmp < percentage) {
    tmp = tmp + vectProportionOfVariance[i]
    i = i + 1
  }
  return(i)
}
nbrPC25 <- numberOfPCs(vectProportionOfVariance, 0.25)
nbrPC50 <- numberOfPCs(vectProportionOfVariance, 0.5)
nbrPC75 <- numberOfPCs(vectProportionOfVariance, 0.75)
nbrPC99 <- numberOfPCs(vectProportionOfVariance, 0.99)
sprintf("Number of principal components needed in order to cover at least 25 percent of variance : %i", nbrPC25)

[1] "Number of principal components needed in order to cover at least 25 percent of variance : 40"

sprintf("Number of principal components needed in order to cover at least 50 percent of variance : %i", nbrPC50)

[1] "Number of principal components needed in order to cover at least 50 percent of variance : 120"

sprintf("Number of principal components needed in order to cover at least 75 percent of variance : %i", nbrPC75)

[1] "Number of principal components needed in order to cover at least 75 percent of variance : 240"

sprintf("Number of principal components needed in order to cover at least 99 percent of variance : %i", nbrPC99)

[1] "Number of principal components needed in order to cover at least 99 percent of variance : 473"

With the test we have established just before we can see that less than \(1/10\) of the data contains \(25%\) of the variance : So it seems like a good idea to use only \(50\) PC to get most of the information within least data. So we might use that number of PC in order to represent genetic markers. The other numbers give a idea of the number needed if we want to cover more information.

4. PCR Principal Components Regression

a)

# pcaNAm2$x are scores : size 494 * 494 matrix. We take all the ligns but only 250 columns
first250PCs <- pcaNAm2$x[, c(1:250)]
# size 494 * 251 => first column become the longitude
long <- c(NAm2$long, first250PCs)
longMatrix <- matrix(data = long, nrow = 494, ncol = 251)
dataFrameLong <- data.frame(longMatrix)
#dataFrameLong
lmlong = lm(formula = X1 ~ . , data = dataFrameLong)
# size 494 * 251 => first column become the latitude
lat <- c(NAm2$lat, first250PCs)
latMatrix <- matrix(data = lat, nrow = 494, ncol = 251)
dataFrameLat <- data.frame(latMatrix)
dataFrameLat

lmlat = lm(formula = X1 ~ . , data = dataFrameLat)
plot(lmlong$fitted.values,lmlat$fitted.values,col="white")
for (i in 1:npop) {
  #print(names[i])
  lines(lmlong$fitted.values[which(NAm2[,3]==names[i])],
        lmlat$fitted.values[which(NAm2[,3]==names[i])],
        type="p",col=colPalette[i],pch=pch[i])
}

legend("bottomleft",legend=names,col=colPalette,lty=-1,pch=pch,cex=.75,ncol=3,lwd=2)

If we compare it to the map from the question 1, we can clearly say that this map represent very optimistically places where population live. It looks quite messy if we just take a look around Panama. However if we examine the estimated places of the Chipewyan, Cree, Ojibwa, Pima and Huiliche populations, we can very distinctly see that the estimation is pretty accurate.

b)

# import the library needed to use rdist.earth
library("fields")

Loading required package: spam
Loading required package: dotCall64
Loading required package: grid
Spam version 2.1-1 (2017-07-02) is loaded.
Type 'help( Spam)' or 'demo( spam)' for a short introduction 
and overview of this package.
Help for individual functions is also obtained by adding the
suffix '.spam' to the function name, e.g. 'help( chol.spam)'.

Attaching package: ‘spam’

The following objects are masked from ‘package:base’:

    backsolve, forwardsolve

# select the theory latitude and longitude of populations
theoryPosition <- matrix(data = c(NAm2$long, NAm2$lat), nrow = 494, ncol = 2)
#theoryPosition
# vector used to stock the mean of distance for each population
stock <- c()
tmpIndex <- 0
for (i in 1:npop) {
  # vector long estimated for population n°i
  tmpLong <- lmlong$fitted.values[which(NAm2[,3]==names[i])]
  
  # vector lat estimated for population n°i
  tmpLat <- lmlat$fitted.values[which(NAm2[,3]==names[i])]
  
  # matrix of the estimated positions
  estimatedPosition <- matrix(data = c(tmpLong, tmpLat), ncol = 2)
  #print(estimatedPosition)
  
  len <- length(estimatedPosition)/2
  
  # by making unique we assume that two estimated coordinates will never be at the same exact position
  diff <- unique(rdist.earth(x1 = theoryPosition[tmpIndex:(tmpIndex+len), ], x2 = estimatedPosition, miles = F, R = NULL))
  
  # stock the mean of distance from the theorical position in kilometers
  stock[i] <- mean(diff)
  
  tmpIndex <- tmpIndex + len
}
plot(c(1:npop), stock, main = "Mean difference in kilometers between theoric and estimated positions",
     xlab = "N° associated to the population", 
     ylab = "Mean of the distance in kilometers")
linearTest <- lm(stock ~ c(1:npop))
abline(linearTest)

Thanks to the previous graph, we can say that the estimated positions are on average approximately \(1100\) kilometers away from the theorical position. It might seem a lot at first sight because it sounds huge if we just have a human sensible approach. Nonetheless it is not really a lot in reality if we size to the scale of the earth. Indeed, we could just think about the fact that the earth’s circonference is equal to \(40075\) kilometers. If we draw a circle around the theorical position for each population, on average, the estimation of the position will be contained inside the circle.

5) PCR and Cross-Validation

a)

Cross-Validation is a technique used in model selection to better estimate the test error of a predictive model. The principle of K-fold cross validation method is to divide randomly the data into K subsets where K-1 subsets are used as the training sets and the remaining subset is used as a validation subset to compute a prediction error. We repeat this procedure K times by considering each subset as the validation subset and the others subsets as the training sets. The procedure is finished when each subset has been used once as the validation subset. The K results from the folds can then be averaged to produce a single estimation.

The advantage of this method is that all observations are used for both training and validation, and each observation is used for validation exactly once.

labels <- rep(1:10, each=10)
set <- sample(labels, 494, replace = TRUE)
set

  [1]  5  6  1  3  1  5  8 10  3  8  3  8  2 10  2  2  2  7  6  3 10  2  1  9  4  9 10  1  3  7  6  4  3  2  7
 [36]  7  7  8  7  1  2  3  3  4 10  7 10  9  8  6  1  3  8  4  7  4  1  1  5  2  8  5  7  4  7  7  7 10  2  2
 [71]  6  3  8  1  7  6 10  2  4  5  3  8  2 10  5  2  3 10 10  7  2  7  9  3  7  3  6  1 10  1  6  8  9  1  6
[106]  1  8  3  5  3  3  7  3  5  2 10  5  9  2  4  6  8  4  5  5  9  4  2  9 10  9  1 10 10 10  6  9  3 10  1
[141]  4  2  2  7  1  9 10  9  8  3  9  6  2  4  5  5  7  6  4  8  8  7  9  9  6 10  6  7 10  1  9  9  5  2  7
[176]  6  5  8  1  2  6  8  8  7  5  5  9  8 10  1  4  2  3  9  2  1  9  2  4 10  6  1  6  4  1  1  8  5  7  9
[211]  7  7  8  5  6  5  4 10  3  6  2  5  4 10  3  6  5 10  6  5  9 10  4  7  5  5  1  9 10  4  2  3  1  4  9
[246]  7 10  4  8  3  1  1  5  5  8  4  6  3  6  4  7  1  8  3  3  3  5  2  4  5  9  6  3  1  4  1  6  8 10  2
[281]  8  1  5  1  2  7 10  3  8  3  2  4  3  4  9  6  1  1  1  2  8  9  4  1  5 10 10  1  5  4  9  2  2  5  9
[316]  8  5  6  9  7  1  1  7  9  8  6  9 10  2 10  3  5  5  6  3  9  4  5  5 10  5 10 10  1  6 10  8  3  4  2
[351]  9  3  1  2  6  1  4  9  2  8  9  4  2  1 10 10  3 10  8  6  4  5  6  7  8  5  8  9  5  8  9 10  4 10  1
[386]  8  9  4 10  2  3  5  1  8  6 10 10 10  2 10  8  2  1  6  3  3  3  8  5  1  2 10  6 10  5  5  8  8  7  6
[421]  7  1  3  4  4  5  3  2  6  8  9  5  1  5  5  8  7  1  3 10  7  3  5  5  2  1  6  1  3 10  9  7  1  1  8
[456]  7  6  9  2  5  6  1  3  7  6  9  2  3  3  8  8  9  5 10  5  6  5  2  2  2  1  1  3  2  1  1  9  1 10  4
[491] 10  9 10  2

b)

naxes <- 4
predictedCoord <- matrix(data = NA, nrow = 494, ncol = 2)
# Regression for the longitude
pcalong <- data.frame(cbind(long=NAm2[,c("long")], pcaNAm2$x[,1:naxes], set))
reglong <- lm(formula = long ~ . , data = subset(pcalong[, 1:5], set != 1))
# Regression for the latitude
pcalat <- data.frame(cbind(lat=NAm2[,c("lat")], pcaNAm2$x[,1:naxes], set))
reglat <- lm(formula = lat ~ . , data = subset(pcalat[, 1:5], set != 1))
# Prediction
# i = 1
predict_long <- predict(reglong, newdata=data.frame(subset(pcalong, set == 1)))
predict_lat <- predict(reglat, newdata=data.frame(subset(pcalat, set == 1)))
idx = 1
for (indice in 1:494) {
  
  if (pcalong[indice, 6] == 1) {
    
    predictedCoord[indice, 1] = predict_long[idx]
    predictedCoord[indice, 2] = predict_lat[idx]
    idx = idx + 1
  }
  
}
# predictedCoord

for(i in 2:10){
    
    reglong <- lm(formula = long ~ . , data = subset(pcalong[, 1:5], set != i))
    reglat <- lm(formula = lat ~ . , data = subset(pcalat[, 1:5], set != i))
    
    predict_long <- predict(reglong, newdata=data.frame(subset(pcalong, set == i)))
    predict_lat <- predict(reglat, newdata=data.frame(subset(pcalat, set == i)))
    
    idx = 1
    for (j in 1:494) {
      
      if (pcalong[j, 6] == i){
      
        predictedCoord[j, 1] <- predict_long[idx]
        predictedCoord[j, 2] <- predict_lat[idx]
        idx = idx + 1
        
      }
      
    }
    
}

# predictedCoord


# Calculation of the prediction error

schtroumph <- rdist.earth(x1 = theoryPosition, x2 = predictedCoord, miles = F, R = NULL)
vect <- schtroumph[1, ]
mean(vect[1:30])

[1] 2289.334

index_vect <- c(1:npop)
indice_main_vect <- 1
index = 1
tmp = NAm2$Pop[1]
for (popName in NAm2$Pop) {
  
  string = popName
  
  if (string != tmp) {
    index_vect[indice_main_vect] <- index
    indice_main_vect <- indice_main_vect + 1
  }
  
  tmp = string
  index = index + 1
}
index_vect[npop] <- 494
main_vect <- c(1:npop)
indice1 <- 1
for (indice in 1:npop) {
  
  indice2 <- index_vect[indice]
  main_vect[indice] <- mean(vect[indice1:indice2])
  indice1 <- indice2
  
}
stock <- main_vect
stock

 [1]  2289.334  3379.760  3541.025  5984.043  5679.188  5844.677  5826.365  7367.050  6671.858  6058.597
[11]  5781.372  6701.283 11418.379  6149.238  6002.236  7543.447  6431.811  5917.866  6884.382  6313.137
[21]  6973.247  5868.110  6697.199  8413.015  9645.011  5494.095  3544.385

plot(c(1:npop), stock, main = "Mean difference in kilometers between theoric and estimated positions",
     xlab = "N° associated to the population", 
     ylab = "Mean of the distance in kilometers")
linearTest <- lm(stock ~ c(1:npop))
abline(linearTest)

c)

for (naxes in 2:440){
  
}